/**
 * 	    \mainpage
 * 		interpolate_regions.cpp
 *
 * 		\section sec_copy Copyright
 *
 * 		Copyright 2011 Stefan Fruhner <stefan@itp.tu-berlin.de>
 *
 *      This program is free software; you can redistribute it and/or modify
 *      it under the terms of the GNU General Public License as published by
 *      the Free Software Foundation; either version 2 of the License, or
 *      (at your option) any later version.
 *
 *      This program is distributed in the hope that it will be useful,
 *      but WITHOUT ANY WARRANTY; without even the implied warranty of
 *      MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *      GNU General Public License for more details.
 *
 *      You should have received a copy of the GNU General Public License
 *      along with this program; if not, write to the Free Software
 *      Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
 *      MA 02110-1301, USA.
 */
#include <mylibs/CardiacMesh.hpp>
#include <mylibs/Profiler.hpp>
#include <mylibs/myinifiles.hpp>
#include <mylibs/mystring.hpp>
#include <mylibs/lists.h>
#include <mylibs/myio.hpp>
#include <algorithm>
#include <ios>
#include <iostream>
#include <fstream>
#include <sys/stat.h> 	// mkdir
#include <omp.h>

using namespace std;
using namespace mylibs;

#ifdef DEBUG
	#define print(x) (cout << __FILE__ << ":" << __FUNCTION__ <<"() (" <<__LINE__ << ") : " << x << endl)
#else
	#define print(x) (cout << x << endl)
#endif

// global variables
string write_basename = "";
mystring input  = mystring();
mystring output = mystring("pacing_nodes/pacing.nodes.vtx");

bool verbose 		 = false;
bool create_list 	 = false;
int  left_ventricle  = 12; // inner of the left ventricle
int  right_ventricle = 13; // inner of the right ventricle
double fraction 	 = 0.3;
int strategy 		 = 0;
int strategy_ref	 = 0;
string method 		 = "pearls";
mystring constraint = mystring("");
bool ignore_z 		 = false;

void save_vtx_file(mystring output, list<int> &nodes, CardiacMesh mesh){
	FILE *f = fopen(output.c_str(), "w");
	fprintf(f,"%d\n", (int) nodes.size());
//	cout << input << endl;
	if (input.file_base().endswith("_i"))	fprintf(f,"intra\n");
	else fprintf(f,"extra\n");
	for (list<int>::iterator i = nodes.begin(); i != nodes.end() ; i++){
		fprintf(f,"%d\n", (int) *i);
	}
	fclose(f);

	output.file_new_ext(".dat");
	mystring gnu(output);
	f = fopen(output.c_str(), "w");
	for (list<int>::iterator i = nodes.begin(); i != nodes.end() ; i++){
		const Point &pt = mesh.point(*i);
		fprintf(f,"%f %f %f\n", pt.x, pt.y, pt.z);
	}
	fclose(f);

	if (not IO::is_executable("gnuplot")){
		cmdline::warning("If you install 'gnuplot', I'll plot the results for you :)");
		return;
	}

	output.file_new_ext(".dat.gnu");
	f = fopen(output.c_str(), "w");
	fprintf(f,"set terminal png nocrop medium size 800,600 \n");
	fprintf(f,"set output '%s'\n", (gnu+"_3D.png").c_str());
	fprintf(f,"splot \"%s\"\n", gnu.c_str());
	fprintf(f,"set output '%s'\n", (gnu+"_2D.png").c_str());
	fprintf(f,"plot \"%s\"\n", gnu.c_str());
	fclose(f);


	mystring cmd("gnuplot " + output);
	if (verbose) cout << " $ " << cmd << endl;
	int result = system(cmd.c_str());
	if (result != 0){
		print("Exit code was " << result);
		std::exit(result);
	}
}

void find_minmax_dist(Point &minp, double &min, double &max, SurfaceMesh &mesh, list<int>& nodes,const Point &com){
	for (list<int>::iterator i = nodes.begin(); i != nodes.end() ; i++){
		Point pt = mesh.point(*i);
		if (ignore_z) pt.z = com.z;
		double d = (pt - com).abs();

		if (d < min) {
			minp = pt;
			min = d;
		}

		if (d > max) max = d;

//		cout << "[" <<pt <<"] " << min << " " << max << endl;
	}

	TELL(min);
	TELL(max);

	return;
}


/**
 * name: applyClippingBox(BoundingBox &bb, SurfaceMesh &mesh, list<int>& nodes)
 *
 * Clear all nodes from the list 'nodes' which are outside of the
 * clipping box. The list contains only indices of points.
 *
 * @param bb	: 	BoundingBox which is to be applied
 * @param mesh	: 	The SurfaceMesh conatining the full node list with
 * 					coordinates
 * @param node	s: 	The list indices of nodes which is to be cleared by
 * 					nodes outside the bounding box. The operation is
 * 					applied in-place.
 * @return void
 **/
void applyClippingBox(BoundingBox &bb, SurfaceMesh &mesh, list<int>& nodes){
	list<int> tmp = nodes;
	nodes.clear();
	for ( list<int>::iterator it = tmp.begin(); it!= tmp.end(); it++){
		if (bb.inside(mesh.point(*it))) nodes.push_back(*it);
	};
}

int main(int argc, char** argv){

	myIniFiles ini(argc, argv, verbose);

	ini.set_info(toString(argv[0]) + "\n\
	(compiled on " + toString(__DATE__) + " " + toString(__TIME__) + ")\n\n\
Generate a vertex list that suffices the 'pearl necklace' pacing \
strategy for 2D slices of the human heart. To account for correct pacing the \
boundaries of the inner part of the ventricles are used as pacing regions. \
The septum is excluded. The output is written as vtx file. Additionally  \
*.dat files are written in order to generate png images using gnuplot. \
\n\
Other methods:\n\
	min_x : select the smallest x-values of a mesh.\n\
	min_y : \n\
	min_z : \n\
");
	ini.register_param("output"	 ,"o", "Basename of files to be written", output.c_str());
	ini.register_param("input"	 ,"i", "The mesh that contains regions and fiber directions.");
	ini.register_param("right"	 ,"r", "ID of the right cavity.", toString(right_ventricle).c_str());
	ini.register_param("left"	 ,"l", "ID of the left cavity.", toString(left_ventricle).c_str());
	ini.register_param("fraction","d", "[0,1] : fraction of the maximum distance that is used to neglect vertices. 1 means no vertices are used, 0 means all vertices are used. default = 0.3", fraction, -1., 1.);
	ini.register_param("strategy","s", "strategy to choose center of septum : 0 = center of mass, 1 = center of closest distance", strategy, 0, 1);
	ini.register_param("strategy_ref","R", "strategy for the reference nodes : 0 = minimum coordinate only, 1 = table |---| ", strategy_ref, 0, 1);
	ini.register_flag("list"	 ,NULL ,"List of possible regions");
	ini.register_param("method", "m", "method to choose pacing sites : ['pearls', 'min_x', 'min_y', 'min_z', 'septum']", method.c_str());
	ini.register_param("constraint", "c", "Constraint for determination of pacing nodes, e.g. z>5 for a minimum z-coordinate = 5 or z<10 for a maxmimum coordinate = 10, both is possible", "");
	ini.register_flag("ignore_z", "z", "ignore z component when method = 'pearls'");
	ini.register_flag("verbose" ,"v", "Additionally messages on stdout.");

	ini.check();

	verbose 		= ini.exists("verbose");
	input			= ini.read("input", input);
	output     		= ini.read("output", output);
	left_ventricle	= ini.read("left", left_ventricle);
	right_ventricle	= ini.read("right", right_ventricle);
	fraction		= ini.read("fraction", fraction);
	strategy		= ini.read("strategy", strategy);
	strategy_ref	= ini.read("strategy_ref", strategy_ref);
	method 			= ini.read("method", method);
	create_list		= ini.exists("list");
	ignore_z		= ini.exists("ignore_z");
	constraint 		= ini.read("constraint", constraint);

	BoundingBox clippingBox;
	if (ini.exists("constraint")){
		cout << "Constraints given : " << endl << constraint << endl;
		clippingBox.fromString(constraint);
	}

//	------------------------- Error checking -----------------------------------
	if (input.empty()) cmdline::exit("-i [--input] : No input mesh ... exit.");

	// if no extension supplied, the assume *.pts
	if (input.file_ext() == mystring()) input.file_new_ext(".pts");
	// check input mesh
	if (not input.file_exists()) {
		mystring err = input + " does not exist.";
		cmdline::exit(err);
	}

	if (strategy > 0) throw myexception("ATTENTION: strategy > 0 disabled !");

//	-------------------- Error checking finished -------------------------------
	mystring base = output.file_base();
	mystring path = output.file_path();
	base.file_new_ext(".pts"); 	// replace extension

	// creating directory
	if (verbose) cmdline::msg("Creating output directory : mkdir("+path+", 0740)");
	IO::mkdir(path.c_str());

	// talk about parameters
	if (verbose){
		cmdline::section("Parameters");
		TELL(input);
		TELL(output);
		TELL(method);
		TELL(left_ventricle);
		TELL(right_ventricle);
		TELL(create_list)
	}

	// load the input mesh
	CardiacMesh mesh(input, verbose);

	if (create_list){
		if (verbose) cmdline::msg(" - Scanning list of regions");
		myList<int> regions = mesh.scan_for_regions();
		if (verbose) cout << "   ";
		for (myList<int>::iterator i = regions.begin(); i != regions.end(); i++){
			cout << *i << " ";
		}
		cout << "\n\n" << endl;

		exit(0);
	}

	// list of accepted nodes
	list<int> pacing_nodes;

	if (method == "pearls"){

		if (verbose) cmdline::section("Determination of boundary nodes");

		Profiler profiler;
		// load the list of nodes at the boundaries
		list<int> left_boundary_nodes 	= mesh.boundary_nodes(left_ventricle);
		list<int> right_boundary_nodes 	= mesh.boundary_nodes(right_ventricle);

		if (ini.exists("constraint")){
			applyClippingBox(clippingBox, mesh, left_boundary_nodes);
			applyClippingBox(clippingBox, mesh, right_boundary_nodes);
		}

		if (verbose){
			cmdline::ok();
			cout << " - "<<left_boundary_nodes.size() << " in the left ventricle" << endl;
			cout << " - "<<right_boundary_nodes.size() << " in the right ventricle" << endl;
			profiler.info();
		}
		profiler.reset();

		if (verbose) cmdline::section("Removing nodes at the septum");
		Point com;	// center of mass

		// find a point in the septum
		if (strategy == 0){
			for (list<int>::iterator i = left_boundary_nodes.begin(); i != left_boundary_nodes.end() ; i++)
				com += mesh.point(*i);
			for (list<int>::iterator i = right_boundary_nodes.begin(); i != right_boundary_nodes.end() ; i++)
				com += mesh.point(*i);

			com /= (left_boundary_nodes.size() + right_boundary_nodes.size());
		} else {
			throw myexception("strategy 1 disabled !");
//			com = (min_l + min_r )/ 2.; // center of shortest distance
		}

// 		determine minimum/maximum distance to the center of mass
		Point min_l = mesh.point(left_boundary_nodes.front());
		Point min_r = mesh.point(right_boundary_nodes.front());
		double minl = 1.E99;
		double minr = 1.E99;
		double maxl = 0.;
		double maxr = 0.;
		#ifdef _OPENMP
			#pragma omp parallel sections
		#endif
		{
	        #ifdef _OPENMP
	        #pragma omp section
	        #endif
	        { // Aufgabe 1
				find_minmax_dist(min_l, minl, maxl, mesh, left_boundary_nodes, com);
			}
			#ifdef _OPENMP
			#pragma omp section
			#endif
			{ // Aufgabe 2
				find_minmax_dist(min_r, minr, maxr, mesh, right_boundary_nodes, com);
			}
		}

		// set overall maxmimum
		double max = (maxr > maxl ) ? maxr : maxl;
		double min = ((min_l - com).abs() < (min_r - com).abs()) ? (min_l - com).abs() : (min_r - com).abs() ;

		if (verbose){
			TELL(com);
			TELL(min);
			TELL(max);
			TELL(min_l);
			TELL(min_r);
			profiler.info();
		}


		// compute allowed distance
		double dist = min + fabs(fraction) * (max - min);

		list<int> LV_nodes;
		list<int> RV_nodes;
		// select nodes for pacing
		#ifdef _OPENMP
			#pragma omp parallel sections
		#endif
		{
	        #ifdef _OPENMP
	        #pragma omp section
	        #endif
	        { // Aufgabe 1
	        for (list<int>::iterator i = left_boundary_nodes.begin(); i != left_boundary_nodes.end() ; i++){
				Point pt = mesh.point(*i);
				if (ignore_z) pt.z = com.z;
				if (fraction > 0.){
					if ((pt - com).abs() >= dist) LV_nodes.push_back(*i);
				} else{
					if ((pt - com).abs() <= dist) LV_nodes.push_back(*i);
				}
			}
	        }
	        #ifdef _OPENMP
	        #pragma omp section
	        #endif
	        { // Aufgabe 2
			for (list<int>::iterator i = right_boundary_nodes.begin(); i != right_boundary_nodes.end() ; i++){
				Point pt = mesh.point(*i);
				if (ignore_z) pt.z = com.z;
				if (fraction > 0.){
					if ((pt - com).abs() >= dist) RV_nodes.push_back(*i);
				}else {
					if ((pt - com).abs() <= dist) RV_nodes.push_back(*i);
				}
			}
			}
		}





		// save RV and LV nodes independently
		mystring outLV = output.file_strip_ext() + "_LV.vtx";
		save_vtx_file(outLV, LV_nodes, mesh);
		mystring outRV = output.file_strip_ext() + "_RV.vtx";
		save_vtx_file(outRV, RV_nodes, mesh);

		// join them in one list
		pacing_nodes.insert(pacing_nodes.end(), LV_nodes.begin(), LV_nodes.end());
		pacing_nodes.insert(pacing_nodes.end(), RV_nodes.begin(), RV_nodes.end());

		cmdline::ok();
	} else if (method == "septum"){

		if (verbose) cmdline::section("Determination of boundary nodes");

		list<int> left_boundary_nodes 	= mesh.boundary_nodes(left_ventricle);
		list<int> right_boundary_nodes 	= mesh.boundary_nodes(right_ventricle);

		if (ini.exists("constraint")){
			applyClippingBox(clippingBox, mesh, left_boundary_nodes);
			applyClippingBox(clippingBox, mesh, right_boundary_nodes);
		}

		if (verbose){
			cout << " - "<<left_boundary_nodes.size() << " in the left ventricle" << endl;
			cout << " - "<<right_boundary_nodes.size() << " in the right ventricle" << endl;
		}

		double min = (mesh.point(left_boundary_nodes.front()) - mesh.point(right_boundary_nodes.front()) ).abs();
		double max = (mesh.point(left_boundary_nodes.front()) - mesh.point(right_boundary_nodes.front())).abs();

		Point min_l = mesh.p[0];
		Point min_r = mesh.p[0];

		for (list<int>::iterator i = left_boundary_nodes.begin(); i != left_boundary_nodes.end() ; i++){
			const Point &pt = mesh.point(*i);

			for (list<int>::iterator j = right_boundary_nodes.begin(); j != right_boundary_nodes.end() ; j++){
				const Point &pt2 = mesh.point(*j);
				double d = (pt - pt2).abs();
				if (d < min) {
					min_l = pt;
					min_r = pt2;
					min = d;
				}

				if (d > max) max = d;
			}
		}

		if (verbose){
			TELL(min);
			TELL(max);
			TELL(min_l);
			TELL(min_r);
		}

		Line line(min_l, min_r, Line::segment);
		Point center = line.point_at(0.5);
		size_t idx = mesh.find_closest_node(center);
		pacing_nodes.push_back(idx);

		cmdline::ok();
	} else if (method == "min_x"){
		Point min = mesh.point(0);
		for (size_t i = 0; i < mesh.points(); i++){
			if (mesh.point(i).x < min.x) min = mesh.point(i);
		}

		for (size_t i = 0; i < mesh.points(); i++){
			if (mesh.point(i).x <= min.x) pacing_nodes.push_back(i);
		}

	}else if (method == "min_y"){

		Point min = mesh.point(0);
		for (size_t i = 0; i < mesh.points(); i++){
			if (mesh.point(i).y < min.y) min = mesh.point(i);
		}

		for (size_t i = 0; i < mesh.points(); i++){
			if (mesh.point(i).y <= min.y) pacing_nodes.push_back(i);
		}

	}else if (method == "min_z"){

		Point min = mesh.point(0);
		for (size_t i = 0; i < mesh.points(); i++){
			if (mesh.point(i).z < min.z) min = mesh.point(i);
		}

		for (size_t i = 0; i < mesh.points(); i++){
			if (mesh.point(i).z <= min.z) pacing_nodes.push_back(i);
		}

	} else throw myexception ("method unknown.");

//	if (output.empty()) output = "pacing_nodes.vtx";
	if (verbose) TELL(output);
	save_vtx_file(output, pacing_nodes, mesh);

	if (verbose) cmdline::section("Determination of reference nodes (for bidomain computations)");

	// find a reference node
	if (strategy_ref == 0){ // minimum coordinate
		if (verbose)cmdline::msg(" - Finding minimum coordinate.");

		// method 1 : search for smallest coordinate in whole mesh
		Profiler profiler;
		size_t minidx = mylibs::algorithm::min_element( mesh.p );
		profiler.info();

		// method 2 : find nodes at the outer boundary and search for the smallest point
//		profiler.reset();
//		list<int> nodes = mesh.boundary_nodes();
//		if (verbose)
//			cmdline::msg(" - Finding minimum coordinate 2/2.");

//		Point min = mesh.point(nodes.front());
//		int minidx = nodes.front();
//		for (list<int>::iterator i = nodes.begin(); i != nodes.end() ; i++){
//			const Point &pt = mesh.point(*i);
//			if (pt < min) {
//				min = pt;
//				minidx = *i;
//			}
//			profiler.info();
//		}

		mystring out(output);
		out = output.file_replace_ext(".vtx", ".reference.vtx");
		list<int> extra;
		extra.push_back(minidx);
		if (verbose) TELL(out);
		save_vtx_file(out, extra, mesh);
		cmdline::ok();
	}else if (strategy_ref == 1){ // strategy == table : 4 grounding nodes
		cmdline::msg(" - Determining 4 grounding nodes");
		Point m;
		Point M;
		mesh.bounding_box(m,M);

		Point p1(m.x, m.y, m.z);
		Point p2(M.x, m.y, m.z);
		Point p3(M.x, M.y, m.z);
		Point p4(m.x, M.y, m.z);

		size_t i1 = mesh.find_closest_node(p1);
		size_t i2 = mesh.find_closest_node(p2);
		size_t i3 = mesh.find_closest_node(p3);
		size_t i4 = mesh.find_closest_node(p4);
		list<int> extra;
		extra.push_back(i1);
		extra.push_back(i2);
		extra.push_back(i3);
		extra.push_back(i4);

		mystring out(output);
		out = out.file_replace_ext(".vtx", ".reference.vtx");
		if (verbose) TELL(out);
		save_vtx_file(out, extra, mesh);
		cmdline::ok();
	}else {
		cmdline::warning("No valid strategy given to determine reference nodes !");
		return -1;
	}

	mystring out(output);
	out = out.file_replace_ext(".vtx", ".ini");
	ini.save(out);

	return 0;
}
